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We have found a new way to express the solutions of the RSM (Reynolds Stress Model) equations that 
allows us to present the turbulent diffusivities for heat, salt and momentum in a way that is considerably 
simpler and thus easier to implement than in previous work. The RSM provides the dimensionless mixing 
efficiencies F * (a stands for heat, salt and momentum). However, to compute the diffusivities, one needs 
additional information, specifically, the dissipation e. Since a dynamic equation for the latter that includes 
the physical processes relevant to the ocean is still not available, one must resort to different sources of 
information outside the RSM to obtain a complete Mixing Scheme usable in OGCMs. 

As for the RSM results, we show that the Ff s are functions of both Ri and R p (Richardson number and 
density ratio representing double diffusion, DD); the F a are different for heat, salt and momentum; in the 
case of heat, the traditional value F u = 0.2 is valid only in the presence of strong shear (when DD is inop- 
erative) while when shear subsides, NATRE data show that F h can be three times as large, a result that we 
reproduce. The salt F s is given in terms of I"/,. The momentum F m has thus far been guessed with differ- 
ent prescriptions while the RSM provides a well defined expression for F m {Ri,R p ). Having tested F h , we 
then test the momentum r m by showing that the turbulent Prandtl number r m /.T| 1 vs. Ri reproduces the 
available data quite well. 

As for the dissipation s, we use different representations, one for the mixed layer (ML), one for the ther- 
mocline and one for the ocean’s bottom. For the ML, we adopt a procedure analogous to the one success- 
fully used in PB (planetary boundary layer) studies; for the thermocline, we employ an expression for the 
variable eN~ 2 from studies of the internal gravity waves spectra which includes a latitude dependence; 
for the ocean bottom, we adopt the enhanced bottom diffusivity expression used by previous authors 
but with a state of the art internal tidal energy formulation and replace the fixed F a = 0.2 with the 
RSM result that brings into the problem the Ri, R p dependence of the T a ; the unresolved bottom drag, 
which has thus far been either ignored or modeled with heuristic relations, is modeled using a formalism 
we previously developed and tested in PBL studies. 

We carried out several tests without an OGCM. Prandtl and flux Richardson numbers vs. Ri. The RSM 
model reproduces both types of data satisfactorily. DD and Mixing efficiency r h (Ri,R p ). The RSM model 
reproduces well the NATRE data. Bimodal £-distribution. NATRE data show that e{Ri < 1) ss 10 s[Ri > 1), 
which our model reproduces. Heat to salt flux ratio. In the Ri > 1 regime, the RSM predictions reproduce 
the data satisfactorily. NATRE mass diffusivity. The z-profile of the mass diffusivity reproduces well the 
measurements at NATRE. The local form of the mixing scheme is algebraic with one cubic equation to 
solve. 
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1. Introduction 

In two previous studies (Canuto et al., 2001, 2002, cited as I and 
II), two vertical mixing schemes for coarse resolution OGCMs 
(ocean general circulation models) were derived and tested. How- 
ever, because of shortcomings in I, II of both physical and structural 
nature, a new mixing scheme became necessary which we present 
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here. By structural we mean that the expressions for the heat, salt 
and momentum diffusivities in 1, 11 were rather cumbersome. By 
physical, we mean the need to include important physical pro- 
cesses that were missing in I, II. 

Concerning the structural issue, we have found a new solution 
of the Reynolds Stress Model, RSM, that yields expressions for the 
diffusivities that are simpler and thus easier to code than the ones 
in II. If we denoted by K a the diffusivities for momentum, heat and 
salt (subscript a), the new solutions of the RSM are: 


Mixed layer : 

2 K 2 

I< a = S x — , 

(la) 

Deep Ocean : 

= r a = hrN) 2 S x 

N 2 

(lb) 


Here, I< is the eddy kinetic energy, e its rate of dissipation, N is the 
Brunt-Vaisala frequency with N 2 = -gp~'p z , t = 2/(e~ 1 is the 
dynamical time scale and S a are dimensionless structure functions 
which are functions of: 


S lx (Ri , Rp,rN) 


(2a) 


where the Richardson number Ri and the density ratio R p (charac- 
terizing double diffusion DD processes) are defined as follows: 



R P 


a. s 8S/dz 
a T 8T /6z 


(2b) 


Here, the variables T, S and U represent the mean potential temper- 
ature, salinity and velocity. The thermal expansion and haline con- 
traction coefficients a TtS = (-p _1 0p/0T,+p _1 6p/0S) may be computed 
using the non-linear UNESCO equation of state and X = (2S,jS,j) 1/2 is 
the mean shear with 2Sy = U !t/ + U,„ where the indices i,j = 1,2,3 and 
a ti = da/dx ,. Relations (la) and (lb) contain two unknown variables, 
the dissipation e and the eddy kinetic energy K: 


which means that to complete the RSM, one must add two more 
relations that provide the variables (3). In engineering flows, these 
two variables are traditionally obtained by solving the so-called I(-e 
model which means two differential equations for those two vari- 
ables. The solution of the K-e model, represented by Eq. (20), would 
close the problem since every variable would now be expressed in 
terms of the large scale fields. Let us analyze how these two vari- 
ables are determined in the present oceanic context. 


1.1. Determination of r 

Since most of the ocean is stably stratified, the vertical extent of 
the eddies is much smaller than the vertical scale of density varia- 
tion (except of course in deep convection places), a local approach 
to the kinetic energy equation, first relation in Eq. (20), is a sensible 
one. Physically, this is equivalent to taking production equal dissi- 
pation, P = e, where P = P s + P b is the total production due to shear 
and buoyancy. Since P = I< m X 2 - I< p N 2 , the derivation is presented 
in Eqs. (22), (23), (54) and (55), use of relations (lb) in P = e trans- 
forms the latter into an algebraic equation for the variable t given 
by Eqs. (40) and (41) the result of which is the function: 

r = T(Ri,R p ) (4) 

Use of (4) in the second of (lb) and in (2a) yields the structure func- 
tions and the mixing efficiencies in terms of the large scale variables: 

S x (Ri,R P ), r a (Ri,R p ) ( 5 ) 

Let us note that the above procedure applies in principle to the 
mixed layer, the thermocline and the ocean bottom. The problem 
is to know how to determine the Richardson number in each region, 


a problem we discuss in Sections 6 and 7.3. When applied to the 
mixed layer, the above determination of the mixing efficiencies is 
physically equivalent to assuming that the external wind directly 
generates oceanic mixing. There is, however, a second possibility, 
namely that the wind first generates surface waves which then be- 
come unstable and break, generating mixing (Craig and Banner, 
1994; Umlauf and Burchard, 2005). To account for such a process, 
one needs the full /(-equation in (20) with a non-zero flux F K of K 
for which one needs a closure. The K-flux F K is a third-order mo- 
ment and, as discussed in Cheng et al. (2005), there is still a great 
deal of uncertainty on how to close such higher-order moments. 
The wave breaking phenomenon is introduced into the problem 
by taking the value of F K at the surface z = 0 equal to the power pro- 
vided by the wave breaking model, as described in the two refer- 
ences just cited. In the present case, local limit P= e, relations (5) 
are still not sufficient to determine the diffusivities given by the 
first relation in (lb) for we require the dissipation e whose determi- 
nation we discuss next. 


1.2. Determination of e 

In principle, one could solve the second of Eq. (20) and obtain 
the dissipation e(Ri,R p ) in analogy with the procedure that lead 
to relations (5). Regrettably, such a procedure is not feasible since 
the equation for e has been problematic since the RSM was first 
employed by Mellor and Yamada (1982). The reason is that, con- 
trary to the /(-equation whose exact form can be derived from tur- 
bulence models, the e-equation has thus far been entirely 
empirically based and a form that includes stable stratification, 
unstable stratification and double diffusion, does not exist in the 
literature. Recently, some progress has been made in deriving an 
e-equation from first principles (Canuto et al., 2010) but only for 
the case of unstable stratification, while most of the ocean is stably 
stratified. For these reasons, we still cannot employ the dynamic 
equation for e and we must rely on a different approach. As for 
the mixed layer, we shall employ the length scheme discussed in 
Section 6, leading us to relations (62)-(64). 1 In the thermocline, 
we borrow from the IGW (internal gravity waves) studies-parame- 
terizations by several authors (Polzin et al., 1995; Polzin, 1996; 
Kunze and Sanford, 1996; Gregg et al., 1996; Toole, 1998) the form 
of e, more precisely, of eN~ 2 , that contains the dependence on lati- 
tude given by Eqs. (65)-(68) which should lead to a sharper tropical 
thermocline. As for the ocean bottom, first we include the enhanced 
bottom diffusivity due to tides, Eq. (70) as suggested by previous 
authors but with the latest representation of the function E(x,y ) (Jay- 
ne, 2009), as well as relation (5) instead of the value r = 0.2 used in 
all previous studies (St. Laurent et al., 2002; Simmons et al., 2004; 
Saenko and Merryfield, 2005); second, the tidal drag given by Eq. 
(72) contains a tidal velocity which thus far has been taken to be a 
constant while we suggest it should be computed consistently with 
the same tidal model that provides the function E(x,y), as we explic- 
itly discuss in the lines after Eq. (72); third, the component of the ti- 
dal field not aligned with the mean velocity cannot be modeled as a 
tidal drag. Since its mean shear is large, it gives rise to a large unre- 
solved shearwith respect to the ocean’s bottom. This process, which 
lowers the local Ri below Ri = 0(1) allowing shear instabilities to en- 
hance the diffusivities, was recognized only in one work by Lee et al. 
(2006) who employed an empirical expression for it. Rather, we 
adopt the knowledge we acquired in dealing with the same problem 
in the PBL (Cheng et al., 2002) which gives rise to relation (73) which 
was tested and assessed in previous work and which was shown to 
work pretty well. 


1 The 1D-GOTM ocean model (Burchard, 2002) has included and solved the e- 

equation in the mixed layer. 
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1.3. Determination of Ri(cr) 

It is part of any RSM to determine whether there is a critical 
Ri(cr) above which mixing vanishes, as it was assumed in the liter- 
ature for many years. The Mellor and Yamada (1982) model pre- 
dicted Ri(cr) = 0.19 which was shown to be so low that the 
resulting mixed layer depths were far too shallow to be acceptable 
(Martin, 1985). In our opinion, the MY result was a motivation for 
the KPP model (which is not based on a turbulence closure) since 
its authors believed that turbulence based models could not give 
better results. Model I yielded an Ri( cr) not 0.19 but 0(1 ) and more 
recently (Canuto et al., 2008a) we showed that there is no Ri( cr) at 
all, as several data of very different nature have now established 
beyond any reasonable doubt. In particular, the new data have 
shown that while the heat flux still decreases toward zero at 
Ri> 1, the momentum flux does not, which means that the surface 
wind stresses are transported deeper than in models with Ri(- 
cr) = 0(1 ). Before using the new mixing scheme in a coarse resolu- 
tion OGCM, and in the spirit of previous schemes such as KPP 
(Large et al., 1994), we carried out a series of tests without an 
OGCM which we briefly describe below. 

(1) In the presence of strong shear, the model predicts that heat 
and salt diffusivities become identical, as expected. 

(2) The model predicts that salt fingers become prevalent at a 
critical density ratio R p « 0.6, in agreement with 
measurements. 

(3) Momentum diffusivity K m (Ri,R p ); most mixing schemes (e.g., 
the KPP model, Large et al., 1994) employ heuristic argu- 
ments. Though lack of direct data does not allow a direct 
assessment of the model prediction of this variable, the pre- 
dicted Prandtl number a t (= ratio of momentum to heat dif- 
fusivities) is shown to reproduce well the measured data vs. 
Ri for the no-DD case, Fig. 3c. Most OGCMs assume a t = 10 
which corresponds to Ri = 0(1 ). 

(4) On the basis of temperature microstructure measurements, 
it was generally assumed that r h = 0.2. Using data from 
NATRE, St. Laurent and Schmitt (1999), cited as SS99, have, 
however, shown that such a value is valid only in regions 
of strong shear and no double diffusion. In the opposite 
regime of weak shear and strong DD, r h can be 3-4 times 
larger. The RSM results yield a r h [Ri,R p ) that fits the data, 
Fig. 4, quite well. 

(5) NATRE data have revealed a bimodal distribution of the 
energy dissipation rate e: in the high e, shear dominated 
Ri < 1 regime, the dissipation is an order of magnitude larger 
than in the low e, salt finger dominated Ri > 1 regime, a fea- 
ture that we reproduce reasonably well, Fig. 5a. 

(6) The heat to salt flux ratio r{Ri,R p ) for Ri > 1 reproduces well 
the values measured at NATRE, as well as laboratory mea- 
surements, Fig. 5b. 

(7) The profile of the mass diffusivity K p at NATRE reproduces 
well the measurements, Fig. 9. 

(8) In the thermocline, in locations where there is no Double 
Diffusion, the RSM model predicts that for Ri(bg) = 0.5 we 
have r h = r s = 0.2, r m = 0.6; while the first two relations 
are as expected, the momentum mixing efficiency turns 
out to be three times as large as those of heat and salt. 

In summary, the complete Mixing Scheme is a combination 
of results from the RSM which lead to a new determination 
of the mixing efficiencies (5) plus prescriptions of how to com- 
pute the dissipation e, the latter being different in different 
parts of the ocean. It is only by combining these two parts that 
one obtains a complete mixing scheme that can be used in an 
OGCM. 


2. Overview of previous and present mixing models 

Ocean general circulation models (OGCMs) solve the dynamic 
equations for the mean (potential) temperature T, salinity S and 
velocity U: 

^(T,S) + U i 0 i (T,S) = -^ j (h i 0,u i s) 

g[j. g (6) 

+ UjdjUj + 2Sij k QjU k = —p 0 6jP — — UiUj 

Here, 0, s, u, are the fluctuating components of the temperature, 
salinity and velocity fields, £1 is the Earth’s rotation, P is the mean 
pressure and Ey k is the totally antisymmetric tensor; overbars de- 
note ensemble averages. To solve Eq. (6), the temperature, salinity 
and momentum fluxes u,0, ujs, ujuj, representing unresolved pro- 
cesses, must be parameterized in terms of the resolved mean vari- 
ables T, S, U. In Eq. (6), the mean velocity field is assumed to be 
incompressible (divergence free) but a treatment of compressible 
flows is available (Canuto, 1997). 

Historically, it was Leonardo da Vinci who, by watching the riv- 
er Arno in Florence, described the water flow as being made of two 
distinct parts, 2 which in modern language are called the mean flow 
and the turbulent, fluctuating component. Several centuries later, 
Reynolds (1895) suggested splitting the total fields into mean and 
fluctuating parts, such as T + 0, S + s, U + u, in what has become 
known as the Reynolds decomposition. The non-linear terms in the 
momentum and temperature (salinity) equations then give rise to 
the terms on the rhs of Eq. (6). Historically, it took a long time to 
realize that Eq. (6) were not the last step of the process. By subtract- 
ing (6) from the equations for the total fields, one obtains the equa- 
tions for the fluctuating fields and from them, one proceeds to derive 
the dynamic equations for the three second-correlations that appear 
in (6). However, such a suggestion was not made until the twenties 
by the Russian mathematician A. Friedmann (the same of the 
expanding universe solution of Einstein’s general relativity equa- 
tions). But, as we shall see in Section 3, even his suggestion was 
not taken up in a concrete form until 1945 (Chou, 1945). Soon after 
0. Reynolds’ proposal, Boussinesq (1877, 1897, cited in Monin and 
Yaglom, 1971, vol. I, Section 3) was the first to suggest heuristic, 
down-gradient type expressions of the form: 

— - 8T SS 0U 

we= ~ Kh dz’ ws= ~ Ks di’ wu = _/Cm 0z (7) 

in which K hfSm represent “turbulent diffusivities”. Several com- 
ments are needed concerning (7). First, even though we have not 
written out the z-dependence explicitly, each function in (7) is com- 
puted at the same z, which means that the model is local. Even 
without knowing the explicit form for the diffusivities (which Bous- 
sinesq did not), it is clear that when large eddies are present, as in 
an unstably stratified, convective region, it is unrealistic to assume 
that the fluxes at a given z are governed only by what occurs in the 
vicinity of z since in reality large eddies span much larger extents so 
large in fact as to be of the same size H of the region, a variable that 
ought to appear in a non-local version of (7), as we show in Eq. (17). 

Stated differently, since by Taylor expansion, to express a non- 
local function one needs an infinite number of derivatives, taking 
only the first of them, as in (7), may not be applicable to convective 
regimes, a topic we shall return to at the end of this section. For the 
time being, however, we assume that locality is an acceptable 
approximation since the majority of the ocean is stably stratified 

2 Observe the motion of the water surface, which resembles that of hair, that has 
two motions: one due to the weight of the shaft, the other to the shape of the curls; 
thus, water has eddying motions, one part of which is due to the principal current, the 
other to the random and reverse motion (translated by Prof. U. Piomelli, University of 
Maryland, private communication, 2008). 
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and the eddies are correspondingly small. This is likely to be the 
reason why the relations (7) have been widely used and are 
amended only when applied to unstably stratified, convective re- 
gimes, as discussed in Section 3. 

The second problem concerns the construction of the diffusivi- 
ties themselves which we have denoted by I< a . Since diffusivities 
have dimensions of (length) 2 time -1 , on dimensional grounds 
alone, one has several relations to choose from: 

K x ~ fV 1 ~ Kx ~ FCV 1 - e 1/3 ^ /3 (8) 

where £ is a typical eddy size, K is the eddy kinetic energy, x = 2/Ce -1 
is the dynamical time scale and e is the rate of dissipation of K. It is 
important to note that the last relation in (8), which follows from 
the preceding one using Kolmogorov’s law K ~ e 2/3 £ 2/3 , was actually 
discovered experimentally by Richardson (1926) 15 years before 
the appearance of the Kolmogorov’s law (Kolmogorov, 1941). How- 
ever, since contrary to the atmospheric related studies of Richard- 
son in which £ represented the separation of two “puffs”, in a 
fully turbulent regime the prescription of l is not straightforward, 
the most physical representation is the third one that involves K, 
e which are calculable quantities for which there exist two dynamic 
equations, Eq. (20). There is a further reason that can be gleaned 
from the definitions of K, e in terms of the spectrum £(k) of the 
kinetic energy: 

K = j E(k) dk, e = 2v J k 2 E(k) dk (9) 

These relations show that /Os integrand peaks at low wavenumbers 
(large scales) while e’s peaks at large wavenumbers (small scales) 
and thus a K-s representation catches both large and small scales. 
It may be useful to recall that even though the second relation in 
(9) contains the kinematic viscosity v, it is known that e is indepen- 
dent of it (Frisch, 1995). Postponing the discussion of how to com- 
pute K-e for a moment, we return to (7) and choose the third 
relation in (8). This gives rise to the two representations (la) and 
(lb). The dimensionless structure functions S x that differentiate heat, 
salt and momentum diffusivities and the mixing efficiencies r x , were 
first introduced in the literature by Mellor and Yarnada (1982) and 
Osborn (1980), respectively. It is quite difficult to guess the struc- 
ture of S x and/or r a with any confidence. The reason is rather 
simple. One must take into account temperature, salinity and veloc- 
ity or more precisely, their gradients, which are usually represented 
by the Richardson number Ri and the density ratio R p defined in 
Eq. (2b). A key task of any mixing scheme is that of constructing 
the structure functions (2a). In the absence of double diffusion, even 
without knowing the exact form of (2a), the general dependence on 
Ri can be guessed at: since shear is a source of mixing, the larger is 
Ri, the smaller must be the diffusivity. It follows that the structure 
functions S x (Ri) must be decreasing functions of Ri. Such general 
argument is at the basis of the Pacanowski and Philander heuristic 
model (1981, PP) in which K s = K h . However, no heuristic structure 
function has yet been proposed for the momentum diffusivity and 
in most OGCMs, K m is treated as a free parameter, e.g., in the GFDL 
model, it is taken to be K m = 1 cm 2 s -1 (Griffies et al., 2005). One 
could in principle improve on that by using available data on the 
turbulent Prandtl number (Webster, 1964; Gerz et al., 1989;Schu- 
mann and Gerz, 1995; Canuto et al., 2008a, and Fig. 3c): 


<h(Ri) = K £ (10a) 

and obtain an Ri-dependent I< m using the PP and/or KPP models for 
Kh with the additional information that at Ri = 0, we have (Canuto 
and Dubovikov, 1996, Eq. (43e)): 

a t (Ri = 0)=^° = 0.72 (10b) 

where Ba and Ko (=1.66) are the Batchelor and Kolmogorov con- 
stants, respectively. 

When double diffusion processes are included, guessing the 
structure functions (2a) as a function of both Ri and R p using only 
heuristic arguments is almost impossible, and the only alternative 
is to adopt the dynamic model known as the Reynolds Stress Mod- 
el, RSM. After the original work of Chou (1945), within the geo- 
physical context the pioneering work was that of Donaldson 
(1973) and Mellor and Yarnada (1982) who derived the structure 
functions: 

Sa(Ri) (11) 

thus opening the way for non-heuristic derivations of such func- 
tions. Since the RSM contains parameters that enter the closure of 
the pressure correlation terms (a detailed discussion can be found 
in several papers, e.g., Cheng et al., 2002), the state of the art of tur- 
bulent modeling at the time of the MY model was such that the 
resulting structure functions (11) decreased rapidly with Ri and 
above a critical Ri(cr) mixing became negligible. Specifically, the 
MY predicted that: 

Ri( cr) = 0.2 (12) 

a value that three years later Martin (1985) showed to yield too 
shallow a mixed layer (ML). The same study also showed that in or- 
der to reproduce the observed much deeper MLs, a value five times 
as large was required: 

Ri(cr) = 0(1) (13) 

It is fair to say that the apparent inability of the original 1982-MY 
model to produce “more mixing" was a key motivation for the KPP 
model (Large et al„ 1994) which is not based on the RSM but on 
an analogy with mixing in the atmospheric boundary layer. 

In 2001, a mixing scheme using the RSM was proposed (I, listed 
in Table 1) which showed that (13) can be derived from the RSM, 
the reason for the difference with (12) being a more complete clo- 
sure model for the pressure correlations and the ability to compute 
several of the constants that were poorly known at the time of the 
MY model but that more recent turbulence modeling allowed to 
compute, as discussed in I. Thus, the primary achievement of I 
was to restore “confidence” in the ability of the RSM to yield re- 
sults in agreement with empirical relations such as (13) by Martin 
(1985). The model, however, had limitations, the most important 
of which are (Table 1): (a) the solutions of the RSM equations were 
rather complex, (b) double diffusion processes were not included, 
(c) mixing due to tides was missing, and (d) a bottom boundary 
layer BBL model was not included. 

In 2002, a second mixing scheme was proposed (II, listed in Ta- 
ble 1) with the goal to include double diffusion processes while the 
other parts of the model were the same as in I. The remaining 


Table 1 

Key features of G1SS mixing schemes. 


Mixing scheme 

RSM 

DD 

K-s 

Ri(cr) 

T: Mixing efficiency 

Latitude dependent IGW 

Tides 

BBL 

I, 2001 

Complex 

No 

Local 

0(1) 

Ri 

No 

No 

No 

II, 2002 

Complex 

Yes 

Local 

R P 

Ri, R p 

No 

No 

No 

III, present 

Simpler 

Yes improved 

Local 

OO 

Ri, R p improved 

Yes 

Yes 

Yes 


74 


V.M. Canuto et al./ Ocean Modelling 34 (2010) 70-91 


shortcomings of II are therefore: (a) the solutions of the RSM equa- 
tions are more complex than in I, (b) no mixing due to tides, and (c) 
no bottom boundary layer. 

In 2008, two new features were found which needed to be 
incorporated into a mixing model: the non-existence of a critical 
Richardson number (Canuto et al., 2008a) and a better DD model 
so as to reproduce the measurements of the heat mixing efficiency 
r h (Ri,R p ) (Canuto et al„ 2008b). Both features are now included in 
the mixing scheme we present here. In Table 1, we summarize the 
key features of the models that have been worked out thus far. 

Since the additional physical features in III naturally made it 
more complex, it was necessary to solve the RSM equations so as 
to obtain a simpler representation of the results than in I, II. Con- 
cerning this point, we need to clarify an important issue. 

The solutions of the RSM provide the structure functions S^Ri.Rp) 
but not the functions I(-e which must be computed separately. This 
means that in the fourth column in Table 1 one could have used 
a non-local model for K-e in any of the three models described thus 
far, which is how Burchard (2002) carried out extensive studies of 
models l-II by adopting the structure functions of those models 
with a non-local model for K-e. 

As already discussed, the fact that the ocean is mostly stably 
stratified, making locality a legitimate approximation, led us to de- 
cide in favor of a local treatment of the K-e equations. It must, 
however, be remarked that it is not clear how poorly local models 
do in an unstably stratified regime such as Deep Convection. For 
that reason, Canuto et al. (2004a) tested the local K-e model with 
the RSM solutions of II in the Labrador Sea and compared the pre- 
dicted mixed layer depths with both observations and predictions 
of KPP and MY-2.5 models in which the equation for K is non-local. 
Comparing the data in Fig. 1 of the paper just cited and the model 
results displayed in its Figs. 2, 3, 9 and 10a, one concludes that, 
while all models predict too deep a ML, mixing scheme II in spite 


of its local nature, performs better than the two non-local models. 
In addition to the non-locality of the K-e equations, there is an 
equally important missing feature, mixed layer mesoscales and 
sub-mesoscales that are known to re-stratify the ML leading to a 
shallower ML, as we discuss in the Conclusions. At this point, it 
is therefore not entirely clear to us how physically relevant is the 
local vs. non-local nature of the vertical mixing model, a feature 
we plan to study in the future. 


3. Structure of the new mixing scheme 

In describing the new parameterization, we follow the items as 
they appear in Table 1 , fourth row, from left to right. We begin with 
the RSM, Reynolds Stress Model which, as already mentioned, has a 
long history whose first application to shear flows appeared in 
1945 (Chou, 1945). For discussions of the RSM, especially in geo- 
physical problems, we suggest the pioneering work of Donaldson 
(1973) and Mellor and Yamada (1982), and the recent reviews by 
Burchard (2002), Cheng et al. (2002) and Umlauf and Burchard 
(2005). The work of Donaldson (1973) is particularly relevant not 
only for its extensive discussion of the closure of higher-order mo- 
ments in terms of lower-order ones, but because it is the first time 
that “four basic principles” were presented in Section 8.4 of Don- 
aldson (1973) which we briefly enumerate: (1) the model must 
be written in covariant or tensor form (so as to be invariant under 
arbitrary transformation of coordinate systems), (2) the model 
must be invariant under a Galilean transformation, (3) the model 
must have the dimensional properties of the term it replaces, 
and (4) the model must satisfy all the conservation properties 
characterizing the variables in question. How these principles 
helped the closure problem is elucidated by several instructive 
examples in the same Section 8.4. 





Fig. 1 . The structure functions S a for momentum, heat, salt and density, see Eqs. (la), (30)-(41) and (55) are plotted vs. Ri for different R p . 
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Fig. 2. Same as Fig. 1 but for the mixing efficiencies r a defined in Eq. (lb). 


3.1. Local and non-local RSM equations 


Though the mean equations (6) require only the fluxes of heat, 
salt and momentum: 


w6 (heat flux), ws (salt flux), u,Uj 
= T,j (Reynolds stresses) (14a) 

the dynamic equations of the variables (14a) involve three more 
correlations (see Appendix A): 

8 2 (temp, variance), s 2 (salin. variance), 

6s (temp. -salin. correlation) (14b) 


For completeness, we present a brief, hopefully illustrative, example 
of how the RSM equations are derived. One begins with the Rey- 
nolds decomposition whereby the full velocity and temperature 
fields are written as the sum of a mean and a fluctuating part,U + u, 
T+0; next, one averages the resulting equations using u = 0 , 0 = 0 
and subtracts the results from the original equations for the total 
fields. The results of this purely algebraic procedure are the equa- 
tions for the fluctuating u ,9 fields that read as follows: 


Du/ 

W 

D0 


6 (UiUj - UiUj) = - ^ - UjUij + CL-rgfi + v 


0Xj 


0X? 


0 2 0 


- + -^ 8 -^ 6 ) = -^, + ^^ 


0X? 


(15) 


plying the second of (15) by 0, one obtains the equation for the tem- 
perature variance and so on. The physical difficulty, known as the 
closure problem, is represented by t he th ird-order moments, 3 TOMs, 
such as P jd. p~uj and ujujuj;, UiUjd, u f 0 2 which physically represent 
the fluxes of Reynolds stresses, heat fluxes and temperature variance 
that must also be “closed”, that is, parameterized in terms of the sec- 
ond-order moments. How this is done was discussed in detail in Can- 
uto (1992) and more recently in Cheng et al. (2002, 2005) and there 
is therefore no need to repeat the discussion here. However, what 
must be stressed is the non-local effects represented by the TOMs. 
The point is that turbulence not only gives rise to non-zero sec- 
ond-order correlations but it also transports them around, that being 
the meaning, for example, of the term ufuJO that represents the flux 
of “heat fluxes”. Similar interpretations apply to the other TOMs. 
When such transport processes are included, one has a non-local 
model since even if the local gradients are zero at some point, imply- 
ing a zero flux and no mixing on the basis of (7), the non-local terms 
ensure the existence of mixing brought about by the fluxes just 
discussed. To give a concrete and simple example, consider the 
equation for the temperature variance_obtained from the second 
of (15). Using the closure k t 00h = -0 2 /x 0 that was derived and 
justified in the papers just cited, one obtains in the ID and stationary 
limits: 


-^01 1 0 —2 

= -T„W0— --!( W0 2 
oz 2 0z 


(16) 


where k t is the thermometric diffusivity and v/k t is the molecular 
Prandtl number («7 for seawater). Since each fluctuating variable 
has zero average, averaging Eq. (15) yields an identity 0 = 0. To ob- 
tain the equations for the second-order moments (14), one proceeds 
as follows: multiply the first of (15) by 0 and the second by u, and 
add the two; the result is the equation for the heat flux u,-0. Multi- 


This shows that where the mean temperature gradient is zero, the 
temperature variance does not vanish due to the flux of 0 2 


3 Along a streamline of an inviscid fluid one has (Euler’s law) jp0u 2 /0f = -dp/dt 
which leads to \pu 2 + p = const, which shows that the pressure is a second-order 
moment. 
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Fig. 3. (a) Dynamical eddy turnover time x = 2/C/e, specifically, G m = (xZ) 2 vs. Ri for different R p solution of Eq. (41); (b) same as in (a) but for the turnover time, specifically, 
G p = (t N) 2 , with the inclusion of double diffusion processes, see Eqs. (26) and (27); (c) turbulent Prandtl number defined in Eqs. (10) with the effect of DD; (d) flux Richardson 
number defined in Eq. (42) with the effects of DD. In panels (c) and (d) the data correspond to the case without DD processes: Kondo et al. (1978, slanting black triangles), 
Bertin et al. (1997, snow-flakes), Strang and Fernando (2001, black circles), Rehmann and Koseff (2004, slanting crosses), Ohya (2001, diamonds), Zilitinkevich et al. (2007, 
2008 LES, triangles), (Stretch et al., 2001, DNS, five-pointed stars). 


represented by the second, non-local, term. Deardorff (1966) 
parameterized non-locality by adding a counter-gradient term y h : 

, ST dw9 2 rr _i _ 

J = wd=-K h — + y h , y h ~z e -^~TH- wj, (17) 

where the “closure” of y h proposed by Holtslag and Moeng (1991) is 
a simplified form of the z-derivative of wO 2 but one that exhibits an 
essential feature, the extent H which represents, as we discussed 
earlier, the fact that when eddies are as large as the “container”, 
the size of H ought to appear in the equations (* represent fiducial 
values). We have also taken t 0 ~ t ~ K/e. It must be stressed that it 
would be unjustified to adopt the same form (17) of the counter- 
gradient also for the salinity and/or momentum fluxes whose form 
requires modeling the corresponding TOMs, an active field of re- 
search, as discussed in a recent work (Cheng et al., 2005). 

The problem of how to solve the RSM equations for the second- 
order correlations (14), including the non-local terms, was studied 
in a recent paper (Canuto et al., 2005) where it was shown how to 
write the fluxes as the sum of local and non-local terms, see for 
example Eq. (9a) of the cited paper. The strategy here is that of first 
solving the local limits of the RSM equations and then adding the 
non-local TOMs terms once a closure form has been chosen. How- 
ever, since the closure of the TOMs is still an active field of re- 
search, it would be premature to adopt any particular closure now. 


S a = S a [(TN) 2 ,R p ,(Tlj 2 ] (18) 

Since in general, N 2 and I 2 appear separately and not as their ratio 
Ri, Eq. (18) does not exhibit the form (5) which entails only two 
large scale variables, Ri and R p . This dissimilarity between (5) needs 
some comments. First, even if we rewrite (18) in the equivalent 
form: 


S a = S a [i?i,R„(TZ) 2 ] (19) 

the function (r Z) 2 still depends on turbulence via the dynamical 
time scale r and therefore at the level of (18) and (19), the problem 
is not “closed". At this point, one has two choices. 

The first choice is to adopt a non-local I<-e model (for a detailed 
discussion of the K-e model, its applications and the e-equation, 
see Pope (2000, Section 10.4) and Burchard (2002)): 


DK 8Fy _ o Ds 8 F £ 
Dt dz ~ ’ Dt + dz 


^(c,P m +c 3 P b -c 2 s) 


( 20 ) 


where c 32 = 1-44, 1.92 and F Ks are the TOMs representing the fluxes 
of K and e: 


Pk = ^wtijtij, F £ = we (21) 

while the buoyancy and shear production terms are defined as 
follows: 


3.2. Ri and R p dependence of the RSM solutions 

The structure functions derived in paper II, Eqs. (13a)-(15), de- 
pend on three variables: 


P b = g(ct T w9 - a S ws) = -goc T — (I( h - K s R p ) 
P m = -(uwU z + vwV z ) = K m I 2 


( 22 ) 

( 23 ) 
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The sign and magnitude of the coefficient c 3 in Eq. (20) are dis- 
cussed in Burchard (2002). In writing (23), we have anticipated 
the fact, to be proven shortly below, that the solutions of the RSM 
are indeed of the form (7). Once a closure for (21) is chosen, the 
solutions of (20) yield K and e and thus t = 2 /C/e. This was the pro- 
cedure used by Burchard (2002) in the 1D-GOTM ocean model in 
which the following closure for the TOMs was adopted: 


one solves Eqs. (A.6)-(A.10) for w 2 . This simple observation 
has allowed us to obtain solutions that are considerably sim- 
pler than those in paper II. 

4. Explicit form of the new mixing scheme 

4.3. Heat and salt diffusivities 


Fk = ~I<m 


e/c 

0Z ’ 


Fe = 




6e 
6 z 


(24) 


The analytic solutions of Eqs. (A.1)-(A.5) yield the following 
form of the dimensionless structure functions: 


To our knowledge, Eq. (20) have not yet been used in 3D global 
OGCMs. The French 3D ocean code OPA employs the first of Eq. 
(20) and a heuristic representation of e in lieu of the second equa- 
tion in (20). 

The second choice is to adopt a stationary, local model for JC: 

P = e, production = dissipation (25) 


and a model for e, as discussed in Section 6. Anticipating that the 
RSM does give rise to diffusivities of the form (1), P = e becomes 
the following algebraic relation: 


P=e: i(T£) 2 S m -i(TN) 2 S p = l (26) 

where: 


Sp 


Si, - S s R f , 
1 -R P 


(27) 


Using Eq. (19), the solution of (26) yields the desired relation (4): 
(T Zf=f(Ri,R p ) (28) 

and thus the final form of (19) is: 

S a = Sa(Ri.Rp) (29) 

which coincides with Eq. (5), thus explaining the conditions of 
validity of the latter. The explicit form of (28) is obtained by solving 
Eqs. (40) and (41). Of course, after determining t we still need a 
model for e which is discussed in Section 6. 


3.3. New strategy for the solution of the RSM 

The RSM equations (5)— ( 11) of paper II are presented in Appen- 
dix A for several reasons: 


(1) to make this paper self-contained, 

(2) to correct misprints in II, specifically, Eqs. (5) and (9), 

(3) to give directly the ID form which is the one being solved 
(using also a simplified notation w"T" — > wO, T" 2 —> 0 2 , 
w"s" — > ws, s" 2 — > s 2 , T"s" — > 6s), 

(4) Eq. (5, II) for the Reynolds stress can be considerably simpli- 
fied by dropping the second term on the rhs of it since the 
coefficient p, is very close to unity, and by rounding off 
the value of p 2 to 1/2. The p, term added much complexity 
to the solution and yet its contribution was quite small. As 
a result, Eq. (A.6) is simpler than Eq. (5, II) and yet it pre- 
serves the key physical ingredients, 

(5) in paper II, we employed a method of symbolic algebra to 
solve Eqs. (5-11, II) simultaneously. The resulting structure 
functions Eqs. (13a-15, II) were algebraic but cumbersome. 
However, inspection of the ID form given in Appendix A 
reveals that this was not an optimal choice since the first five 
equations do not depend on shear which appears only in Eq. 
(A.6). Thus, one can separate the problem into two parts: 
first, one solves Eqs. (A.1)-(A.5) analytically (without the 
need of symbolic algebra methods) and in a second step, 


Sh,s — A/ IS 


w 2 


(30) 


where: 

A h = 7t 4 [1 + px + 7t 4 7t 2 x(l -r -1 )] -1 , A s = A^rR,,)^ (31) 

Following standard notation, we denote by r the heat-to-salt flux 
ratio given by the following relations: 

r = ctrw0 = j_Kh K h = 7t 4 1 + qx 
~ a s ws R p K s ’ I< s 7t i 1 +px 

where the dimensionless variables x, p and q are defined as follows: 

x = (tN) 2 (1 -R p ) \ p = 7i 4 7r 5 - 7i 4 7r 2 (l +R P ), 
q = 71, 7t 2 (1 + R p ) - 7t, n 3 Rp (33) 

The tiffs are the dissipation - relaxation time scales defined in Eq. 
(12) of II made dimensionless by using the dynamical time scale 
t = 2 K/e. As one can observe, the fact that we have not yet used 
Eq. (A.6) for the Reynolds stresses is manifest in the still unknown 
w 2 /K term in (30) which we determine next. 

4.2. Momentum diffusivity 


Consider the Reynolds stress equations (A.6)-(A.10). In the sta- 
tionary limit, one obtains a set of linear algebraic equations in the 
variable by which can be solved. The structure functions for the 
case of momentum have the following form: 



Am 


A m 1 
A m 2 


where: 


71ml = 
Am2 = 


4 

5 


7l 4 - 71, + ( 7t, - — ](1 -r ') 


xA h 


10+ (7I 4 - 7I,R,,)x + — (xZ) 


(34) 


(35) 

(36) 


4.3. The ratio w 2 /K 


The general form of the ratio w 2 /I< is given by: 


w 2 

T 


1+ 7^ X + l W TZ)2 


, X = (1 - r _1 )xAh 


(37) 


It is important to note that, contrary to what has been done in many 
ocean models, it is no longer necessary to guess the momentum dif- 
fusivity, as discussed in Section 1, since the model provides K m as it 
provides heat and salt diffusivities. In conclusion, the above formu- 
lation shows that all the variables exhibit the dependence on the 
three functions: 


(zN) 2 ,(TE) 2 ,R p ^Ri,Rp,(Tl) 2 


(38) 
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4.4. Dynamical time scale x 


4.5. Overview 


If one solves Eq. (20) for K and e, x = 2 K/e is automatically given 
as a function of Ri and R p . If, on the other hand, one uses the local 
model represented by the assumption P = e, simplifications of the 
above equations are possible. Expressions (31) for A hs remain the 
same while expressions (35)-(37) simplify to: 



1 w 2 

2l 



(39) 


Next, using the notation by Mellor and Yamada (1982): 


The mixing model is now complete since momentum, heat and 
salt diffusivities have been expressed in terms of the resolved fields 
represented by two large scale variables Ri and R p , and Eq. (6) can 
therefore be solved. 

5. Tests of the mixing model without an OGCM 

Before using the above mixing model in an OGCM, we believe it 
is important to assess its validity and predictions without using an 
OGCM. In what follows, we present the tests we have carried out. 


G m = (TZ) 2 


(40) 5 . 1 . Test: strong shear 


Eq. (26) becomes a cubic equation for G m in terms of Ri and R p : 
c 3^m + C 2^m + CiGm + 1 = 0 

C3 = A^Ri A 2 R 1 , C2 = A^Ri + A4RL Ci ^AsRi + As 

The functions A k ’s are given in Appendix B, Eqs. (B.12). Once the 
function G m (Ri,R p ) is known, one can construct all the relevant func- 
tions the most prominent of which, the structure functions, are pre- 
sented in Fig. 1, the mixing efficiencies in Fig. 2 and the time scales 
in Fig. 3a and b. For example, Fig. 3a exhibits several new features; 
in mixing model II with a finite Ri(cr), the function G m became infi- 
nitely large at Ri(cr) = 0(l) corresponding to the vanishing of the 
eddy kinetic energy since x ~ £K~^ 2 , which is the way Ri( cr) was de- 
fined in that scheme. An alternative interpretation is that at Ri( cr), 
the eddy lifetime becomes very large indicating a tendency toward 
laminarity, that is, in the absence of the breakups of the linear struc- 
tures by the non-linear interactions, the eddy life times x -» 00. In 
the present mixing scheme with no Ri( cr), in the case without DD 
processes R p = 0, G m still increases as Ri increases and turbulence 
decreases but it no longer diverges at any Ri reaching instead a fi- 
nite asymptotic value. However, in the presence of DD processes 
and at sufficiently large Ri corresponding to a vanishing shear (a 
source of mixing), the DD itself becomes a source of mixing which 
leads to a decrease of the eddy life time and G m decreases corre- 
spondingly. On the other hand, the results also show that as long 
as shear is strong (small Ri), DD has no effect being overpowered 
by the stronger action of shear and thus all R p give the same result 
as R p = 0. It is known from laboratory data (Linden, 1971) that 
strong shear disrupts salt finger formation. As the data presented 
in Figs. 9 and 10 of Canuto et al. (2008a) show, the lack of an Ri{ cr) 
is more evident in the momentum than in the heat flux which be- 
comes very small at Ri >0(1) and that is why the function G p = 
(xN) 2 in Fig. 3b still grows with Ri when DD processes are not pres- 
ent ( R p = 0), and why the presence of DD processes softens the 
growth but does not have nearly as dramatic an effect as in Fig. 3a. 

In Fig. 3c we present the turbulent Prandtl number, Eq. (10), vs. 
Ri. We have superimposed data for the no-DD case to show that the 
model reproduces them satisfactorily. Finally, in Fig. 3d we plot the 
flux Richardson number derived from Eqs. (22), (23), (27) and (32): 

P = K m Z 2 (l-R f ), R f = Ri^ = Ri^\^ (42) 

*'m I'm 1 ftp 

For the R p = 0 no-DD case, the available data are well reproduced 
since (1 - r _1 )(l - R p )~’ is unity in this case. DD processes affect 
Rf quite significantly: in the strongest DD case considered here, 
Rp = 0.8, Rf becomes negative quite early. The physical interpreta- 
tion of Rf< 0 is that the buoyancy flux, instead of acting like a sink 
as in the absence of DD processes, becomes a source of mixing due 
to salt fingers instabilities and contributes positively to the total 
production P. 


Since x defined in Eq. (33) represents the eddy turnover time t, 
a strong turbulent regime (Ri <c 1 ) corresponds to small x’s and a 
small x, in which case the last relation in Eq. (32), together with 
Eq. (A. 11), yields: 

Kh = I< s (43) 

which is a reassuring result since when mixing is strong, such as in 
the ocean’s wind driven mixed layer, there is no difference between 
salt and heat diffusivities, as it was proven in laboratory experi- 
ments (Linden, 1971). Double Diffusion processes can only operate 
when shear has subsided, which occurs below the ML. 

5.2. Test: weak shear 

This case corresponds to neglecting shear in (26). Using (32) and 
(33), Eq. (26) acquires the form: 

l(TN) 2 S h =(l-R p )(r- 1 -l)-’ (44) 

Using the representation (lb), that is: 

Ry. = r % * , r a = hxN) 2 s a (45) 

iV ^ 

and combining Eqs. (44) and (45), the mixing efficiency for the tem- 
perature field is given by: 

^(l-i^Xr- 1 -!)-’ (46) 

In the case of salt fingers, measured data give r rs 0.6-0.7, R p r= 0.6- 
0.7 (Kunze, 2003; Schmitt, 2003) and thus the model predicts that 
r h has the value: 

r h = 0.6— 0.7 (47) 

in agreement with the last panel in Fig. 4. We also note that (47) is 
more than 3 times larger than the canonical 0.2 with no double dif- 
fusion (Osborn, 1980, Eq. (10)). 

5.3. Test: onset of salt fingers at R p (cr) 

Next, we assess the model ability to predict the value of R p (c r) 
that characterizes the onset of salt fingers (SF) and diffusive con- 
vection (DC). It is important to recall that linear analysis predicts 
that SF occur in the regime (Schmitt, 1994): 

10~ 2 ^ R p ^ 1 (48) 

Kt 

where the lhs is the ratio of the salt to heat kinematic molecular dif- 
fusivities. On the other hand, Schmitt and Evans (1978) showed that 
in the ocean SF become strongly established when: 

R p > R p ( cr) ~ 0.6 (49) 
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Fig. 4. The heat mixing efficiency r h (Ri,R p ) defined in Eq. (lb). The model results (dashes and full lines) are superimposed on the data from NATRE-TOPO from St. Laurent and 
Schmitt (1999, Fig. 9). 


which is quite different than (48). First, we comment on the fact 
that the present model is sufficiently general to encompass (48) 
in the appropriate limit and then show that it does yield the correct 
result (49). The present RSM formalism is valid for arbitrary dissipa- 
tion-relaxation time scales that appear in the last terms in Eqs. 
(A1)-(A5) and (A9) and (A10) and which, when written in units of 
the dynamical time scale t, are denoted by the n’s which in general 
depend on both the kinematic molecular heat, salt and momentum 
diffusivities, a well as on Ri and R p . Relations (A.ll) used here cor- 
respond to relatively large Reynolds numbers Re. In the limit of 
small Re, the form of the n’s was given by Zeman and Lumley 
(1982) and when used in (32), it yields (48). In spite of several at- 
tempts, we have not yet been able to find a general expression for 
the n’s valid for all Re. To show that the present large Re model 
yields a value corresponding to (49), we consider Fig. 2b which rep- 
resents the heat mixing efficiency r h vs. Ri for different R p . The 
interesting feature is that as one begins with R p = 0 and increases 
its value, there is an uppermost curve past which a further increase 
in R p corresponds to lower values of r h . The R p value corresponding 
to the maximum r h , which we shall call R p { cr), can be read from the 
curves to be around: 

R p {c r) ~ 0.6 (50) 

which reproduces (49) corresponding to the onset of SF. 




2 e (0T/8z) 2 


Let us note that (51 ) corresponds to the stationary limit of Eq. (A.3) 
where / = 20 2 t ( ) 1 , with x 0 = rn 5 . SS99 results, shown in Fig. 4, exhi- 
bit new and interesting features, the most prominent of which is the 
fact that the canonical value r h = 0.2 which has been used for years, 
is valid only in the presence of strong shear when double diffusion 
(DD) processes cannot operate. However, when shear subsides and 
DD become active, the mixing efficiency becomes 3-4 times as large 
(Fig. 4f). 

The challenge for any mixing scheme is to reproduce the data of 
Fig. 4. We begin by showing that Eq. (4) of SS99 (we recall that 
R P (SS99) = /?“’): 


r _ Rf 1 Rp 
h ~ 1 - Rf 1 - r-i 


(52) 


is identical to our Eq. (lb) under P = e. First, the flux Richardson 
number, including double diffusion processes, is defined as: 


R f J^Ri = 


r P 

i + r p 


(53) 


The buoyancy diffusivity K p follows from (22) rewritten as: 


5.4. Test: mixing efficiency r h 

Using NATRE and TOPO data to estimate % (rate of dissipation of 
the temperature variance) and e, SS99 plotted the heat mixing effi- 
ciency r h (Oakey, 1985) as a function of Ri and R p : 


P b = -g«T fiz (K h - K S R P ) = - K P N 2 (54) 

where K p is the mass diffusivity given by: 

K p = /<h(l - r _1 )(l - Rp) ’ 


( 55 ) 
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This allows us to write the total production P in the compact form: 

P = I< m Z 2 - K P N 2 (56) 

Eliminating I< m between (53) and (56), the form of r h given in (52) 
becomes: 

N 2 1 , 

r h = T I< h = -(TE) 2 RiS h (57) 

which is Eq. (lb). The model predictions based on Eq. (57) are plot- 
ted in Fig. 4. The model reproduces the major features of the SS99 
data. 

To put the model results of Fig. 4 in perspective, we point out 
that to the best of our knowledge, no mixing model has thus far 
reproduced these data. In fact, all treatments we have seen employ 
a heat mixing efficiency of 0.2 irrespectively of whether there are 
DD processes or not. It is perhaps more accurate to say that they 
neglect DD and in so doing, they underestimate the true mixing 
which, as demonstrated in Fig. 4f, can be up to three times as large 
as the No-DD case. 

5.5. Test: low and high e-modes 

SS99 analysis of the NATRE data further revealed a bimodal dis- 
tribution of the kinetic energy dissipation rate e. The first regime, 
called the high e-mode, is characterized by a dissipation rate of 
the order of 10 -9 W/kg, while the second regime, called the low 
e-mode, is characterized by values 10 times smaller than those of 
the first mode, e = 0.1 x 1CT 9 W/kg. The first mode was identified 
with an Ri < 1 turbulent regime while the latter was identified with 
an Ri> 1, salt finger dominated regime. The challenge posed by 
these data to any mixing model is not trivial. The reason is that 
the latter usually do not include dynamic equations for the dissipa- 
tion variables e, x which in SS99 were taken from the data them- 
selves. An heuristic equation for e exists, second of Eq. (20) but 
lacks double diffusion, thus making the second of (20) too incom- 
plete for the case under study. We therefore suggest the following 
alternative. Using Eq. (51), we derive the following relation: 

S = ~ R ^ (1CT 9 W/kg) (58) 

4 N 2 r h 

where we have used the following notation: N 2 t = N 2 /Nj, ATRE , 
Nnatre = 1 -7 x 10- 5 s- 2 , C = x 9 oc 2 , * 9 = Z/1(T 9 K 2 s- 1 and a 3 = a T / 
(3 x 10 4 K 1 ). Using the mixing model result shown in Fig. 4 for 
the mixing efficiency T h , in Fig. 5a we plot relation (58) for the 


NATRE case corresponding to C = 1 and N 2 = 1. From the figure 
we deduce that: 

e(Ri = 0.05) = 10e(Ri =10) (59) 

a result in general agreement with the SS99 finding. It must be 
noted, however, that contrary to the case of the mixing efficiency 
r h (Ri,R p ) for which the mixing model provided the full function 
r h (Ri,R p ) which we compared directly with the data, in relation 
(58) the uncertainties still present in modeling / are such that 
the mixing model is unable to provide the full function £(Ri,R p ). 
To arrive at the results presented in Fig. 5a, we borrowed the func- 
tion x from the SS99 data. Regrettably, at this stage of model devel- 
opment, we cannot do any better. 

5.6. Test: heat to salt flux ratio r(Ri,p) 

In the Ri » 1 case, the heat to salt flux ratio is given by Eq. (46) 
which we rewrite as: 


which depends on the mixing efficiency r h (Ri,R p ) which we have 
already assessed in Fig. 4. Eq. (60) is plotted in Fig. 5b together with 
the data of Fig. 10a of SS99. The SS99 data, represented by the blue 
line and the gray area representing the errors, are satisfactorily 
reproduced by the model. 

5.7. Test: momentum diffusivity and turbulent Prandtl number 

The turbulent Prandtl number, given by Eqs. (10), (1) and (30)— 
(37), is shown in Fig. 3c. In most OGCMs, the momentum diffusiv- 
ity K m is treated as a free parameter with a value frequently taken 
to be 10 times larger than l< h , that is, 1 cm 2 s \ which means a tur- 
bulent Prandtl number of 10 which corresponds to an Ri = 2-4, as 
shown in Fig. 3c, which is larger than the value corresponding to 
the internal gravity waves field, as discussed in Section 7.2. 

5.8. Previous models 

In paper II, we reviewed some previous models and here we 
need to update the discussion. We begin with the work of Smyth 
and Kimura (2007) who employed linear stability analysis to study 
DD under the influence of shear. When they compared the model 
results for r h with the data of Fig. 4, the predicted dependence 
on R p was the opposite to that of the data. Inoue et al. (2007) 


SF 




Fig. 5. (a) Plot of the kinetic energy dissipation e in units of 1 0 9 W/kg given by Eq. (58) vs. Ri for R,, = 0.5; the two values at Ri = 0.05 and 10 yield a ratio of 10 which is in 
accord with the SS99 finding of a bimodal distribution of the kinetic energy dissipation; (b) plot of the heat to salt flux ratio r(Ri,R,,) given by Eq. (60), for vanishing shear. The 
NATRE data of SS99 are indicated by the blue line with the gray area indicating the error bars. The model reproduces the data satisfactorily. (For interpretation of the 
references in color in this figure legend, the reader is referred to the web version of this article.) 
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employed a model similar in spirit to the partition first suggested 
by Walsh and Ruddick (2000) which reads: 

K a (Ri,Rp) = K a (Ri > ),R p ) +K a (Ri < 1) (61) 

DD Turb 

with the understanding that the turbulent part no longer depends 
on the density ratio R p . Since the ratio K(turb)//((DD) is not given 
by the WR model, it was treated as a free parameter. Inoue et al. 
(2007) defined the crossing point when the Reynolds number Re = - 
e(vN 2 ) -1 = 20; they further employed a heuristic expression for the 
salt diffusivity K s in the SF regime suggested by Zhang and Huang 
(1998) while for I< h they employed the first of (32) with a constant 
r = 0.71. For the DC regime, they employed a model forKj, and r sug- 
gested by Kelley (1990). However, as none of their relations con- 
tains Ri, it seems unlikely that they can reproduce the data in Fig. 4. 

As for coupled global oceanic-atmospheric codes, the GFDL code 
(Griffies et al., 2005, see their Eqs. (2)-(4)) accounts for SF but not 


DC and employs laboratory data to model SF. However, since there 
is no Ri dependence in such DD model, the resulting heat mixing 
efficiency may be underestimated when SF processes are strong. 

We complete the discussion about DD with some brief remarks 
about their oceanic importance (Ruddick and Gargett, 2003). WR 
noted that at NATRE (Ledwell et al., 1993, 1998) the diapycnal mix- 
ing of heat, salt and tracer is dominated by turbulence but en- 
hanced by salt fingers, and Kelley (2001) noted that at NATRE up 
to half of the diffusion (of an injected tracer) might have been 
transported by DD (St. Laurent and Schmitt, 1999; Kelley, 2001). 
Furthermore, from the maps of the oceanic sites susceptible to 
DD process presented in Figs. 6, 7 (kindly provided to us by Dr. 
D.E. Kelley), one observes that the likelihood of SF (salt fingers) 
processes is higher in the Atlantic (the location of NATRE) than 
in most of the Pacific and that DC (diffusive convection) may play 
a significant role in the Arctic and Southern oceans, a point dis- 
cussed with extensive references by Kelley et al. (2003, Sec- 
tion 2.3.2) who concluded that DC “could be of major importance 


SF 



i 1 1 1 1 1 1 1 r 

Fig. 6. Salt Fingers. Ocean regions susceptible to SF; R~ J intervals are: 1-1.5 (red), 1.5-2 (light brown) and 2-3 (yellow). The reddest color has the lowest which is most 
favorable to SF. Courtesy of D.E. Kelley. 



Fig. 7. Diffusive convection. Ocean regions susceptible to DC; R p intervals are: 1-3 (red), 3-5 (light brown) and 5-10 (yellow). The reddest color has the lowest R p which is 
most favorable to DC. Courtesy of D.E. Kelley. 
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to the properties of the global ocean”. In general, DC is more likely 
in high-latitude precipitation zones (Schmitt, 1994) and Muench 
et al. (1990) also found it in Antarctica over much of the Weddell 
Sea. Overall, in the circumpolar current, both SF and DC may be 
quite important. In conclusion, the results of the present model 
are closer to the data than those of any previous models. 

6. Modeling dissipation 

As previously discussed, Eq. (20) for the dissipation e has never 
been derived from first principles using two-point closures (see, 
however, Canuto et al. (2010)), it contains adjustable parameters 
whose sign is still in dispute, it does not contain double diffusion 
and it is not clear how to extend it to include internal gravity 
waves. Under such circumstances, the best one can do is to employ 
heuristic models, as we now discuss. 

6.1. Mixed layer 

We follow the methodology discussed in paper I, Section 1 1 and 
paper II, Section 9a and model e as follows (Mellor and Yamada, 
1982): 

£ ML = i/fS 3 , // = ri 0 (TZy 3 (62) 

Since the mixing length £ and the shear I are the natural variables, 
the combination £ 2 ! 3 follows; the second relation in (62) comes 
from using the following relations e = K 3I2 IA, t = 2/C/e, A = 8~' 1 
2 B if, where the numerical coefficient t; 0 = B 2 stems from the rela- 
tion A = 8 where B ^ = G 3 J 4 (Ri = 0) = 21.6 as discussed in 

Cheng et al. (2002). As for the mixing length, we employ the 
relation: 

V = (KZ)-'+V (63) 

which follows from Blackadar (1962) who suggested that the mix- 
ing length £ be taken as half of the harmonic mean of kz and 
£ 0 = Q.\1H, H being the depth of the mixed layer and k = 0A the 
von Karman constant. The z-dependence in (63) is such that for 
small z’s, one recovers the law of the wall £ ~ kz, whereas for larger 
z’s, £ becomes a constant fraction of H, as indicated by LES (Moeng 
and Sullivan, 1994). Following previous authors, e.g., Large et al. 
(1997), H is where the potential density differs from the surface va- 
lue by |cr(H) - cr(0)| > 3 x 10~ 5 g citT 3 . However, Zilitinkevich et al. 
(2007) found that to match boundary layer data the length scale 
had to be reduced for large values of the flux Richardson number. 
For the purposes of a model which includes salt and heat contribu- 
tions to stratification, we use the flux Richardson number defined in 
(53). Since the ratio K p II< m depends on Ri, R p , so does Rf. As Fig. 8a 



shows, at each R p , as Ri increases towards infinity, Rf approaches a 
finite limit R y, which is still a function of R p . Generalizing the for- 
mula of Zilitinkevich et al. (2007), we then write: 

The factor introduced by Zilitinkevich in the above length scale 
causes the mixed layer contribution to the diffusivity to fall off 
quickly below the mixed layer so that we may use it in the region 
below as well to allow a continuous transition. In the mixed layer, 
Ri is computed using the resolved large scale fields. 

6.2. Thermocline 

In this region, we have two main contributors, 1GW (internal 
gravity waves) and double diffusion which we now discuss. We be- 
gin by generalizing relation (lb) in the following way: 

K^r^R^WN) (65a) 

with: 

L{0, N) = [fArcosh(N//)] [f 30 Arcosh(7V 0 //3o)] 1 (65b) 

where f 30 means / computed at 30°, N 0 = 5.24 x 1CT 3 s _1 and e t h is 
the dissipation in the thermocline. The function L(0,N) accounts 
for the measured latitudinal dependence of the IGW spectra (Gregg 
et al., 2003) which affects all diffusivities. The effect of the latitude 
dependence on the sharpness of the tropical thermocline was stud- 
ied in Canuto et al. (2004b). As for e th , the best procedure would be 
to identify it with relation (58) which, formally, is quite general. 
This is, however, not a feasible procedure because we lack a model 
for the dissipation x(Ri,R p ) of general validity while we have only a 
few values measured at selected locations, e.g., NATRE. Use of (58) 
everywhere would therefore be unjustified. 

Next, consider the contribution of IGW to £ igw which we quan- 
tify using the Gregg-Henyey-Polzin model (Polzin et al., 1995; Pol- 
zin, 1996; Kunze and Sanford, 1996; Gregg et al., 1996; Toole, 
1998) which gives: 

^ = 0.288A (cgs units) (66) 

where the dimensionless constant A accounts primarily for devia- 
tions from the Garrett-Munk background internal gravity wave 
spectrum and it varies at most by a factor of 2. If we employ again 
the NATRE data with A = 1, N 2 = N^ ATRE , we obtain from (66): 

£igw = 0.5 (!CT 9 W/kg) (67) 



Fig. 8. (a) Same as in Fig. 3d with Rf re-plotted on a linear scale to exhibit negative values; (b) the values of R p corresponding to r p = 0 derived from Fig. 2d. The black dot 
indicates the value of R p = 0.61 past which SF dominate over shear. The corresponding Ri is discussed in the text. 
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which is in the middle of Fig. 5a. For example, with an efficiency of 
r h = 0.25, Eq. (66) yields a diffusivity of 0.07 cm 2 s _1 which is in 
accordance with the NATRE measurements (Ledwell et al., 1993). 
At present, lacking a more complete model, we shall use relations 
(65) with: 

^ = 0.288 (cgs units) (68) 

The last problem concerns the Ri in (65). It cannot be identified with 
the large scale Ri for it would yield practically zero diffusivity. It 
must be a much lower value, which we call Ri (back), which is con- 
tributed mostly by the shear generated by the internal waves which 
is not resolved by the OGCMs and which must therefore be mod- 
eled. To identify Ri (back), we suggest the following procedure. Con- 
sider the plot r p vs. Ri (Fig. 2d): we observe that r p becomes 
negative at different Ri for different R p . The physical meaning of 
the transition is as follows. While shear produces r p > 0, DD pro- 
duces r p < 0: thus, we identify the change from r p > 0 to r p < 0 
as the transition to a DD regime. In addition, since Schmitt and 
Evans (1978) and Zhang and Huang (1998) showed that SF become 
prevalent only at/or past R p sw 0.6, this value is plotted in Fig. 8b 
(Ri - R p points corresponding to r p = 0) as a horizontal line. The 
corresponding Ri ss 0.5 is taken to be the value of Ri (back). 

It must be stressed that we are not suggesting that the mea- 
sured values of Ri below the mixed layer be identified with Ri 
(back). Instead, we view Ri (back) as an effective Ri at which the 
diffusivity approximates the average of the diffusivities over a re- 
gion of space and time containing points with a wide range of Ri. 
We take the point of view that the heat and salt diffusivities pro- 
duced by SF, IGW shear mixing, and the interaction of the two, 
have spatial and temporal scales larger than those of the two pro- 
cesses separately. In building models for coarse OGCMs that do not 
resolve IGW or patches of SF, we can only attempt to model these 
large scale diffusivities. While the offline results in Fig. 4 show that 


our mixing model can reproduce the results of local measure- 
ments, they also illustrate the difficulty of translating such success 
into an OGCM parameterization. In fact, even after restricting to a 
SF favorable R p and removing 75% of the data with lower dissipa- 
tion, SS99 data in Fig. 4 still show a wide range of Ri, that is, for 
a fixedRp, the measured Ri may vary from less than 0.25 to greater 
than 5. A single OGCM point represents a range of conditions 
including those where wave breaking produces strong shears and 
small Ri, as well as quiescent regions where no wave-breaking is 
occurring and Ri is large. With Ri (back), we attempt to represent 
the effects of this whole range of Ri' s that the OGCM does not 
resolve. Although most of the data points in Fig. 4 for R p = 0.6 have 
Ri > 0.5, it must be remembered that the lowest Ri entails large 
diffusivities and thus carry greatest weight in the average 
diffusivity. 

However, there remains the question of the dependence of Ri 
(back) on R p . In paper 11, Ri (back) was taken to be a fraction of 
Ri(cr), the latter being traditionally defined as the value past which 
turbulent mixing vanishes. On the other hand, since in the present 
more realistic model there is no longer an Ri( cr), the approach in II 
is no longer viable. We have examined several alternatives to find 
the R p dependence of Ri (back) but have not yet obtained a credible 
result. While the search continues, we decided to examine the sim- 
plest case of taking Ri (back) constant. Since Ri (back) is introduced 
to produce average diffusivities for coarse resolution OGCMs, it 
must be tested against averaged data. In Fig. 9, we compare the 
mass diffusivity I< p from the model with Ri (back) = 0.5 with the 
SS99 data which are averages over many measurements (the 90 
meter point was excluded because there may be contamination 
from the boundary layer and the error bar on the measurement 
is quite large). 

Since the model vs. data are in reasonable agreement, we feel 
that while we keep on searching to improve it, the relation Ri 
(back) = 0.5 is a tolerable, provisional approximation. 
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Fig. 9. NATRE: the crosses represent the mass diffusivity K p defined in Eq. (55) as a function of depth at the location of NATRE without use of an OGCM. The heat and salt 
diffusivities in Eqs. (55) and (32) are given by the model as a functions of both Ri and R p , together with Eq. (68). The values of R p are taken from the St. Laurent and Schmitt 
(1999) data while Ri is taken to be 0.5, as discussed in the text. The error bars of the model results reflect the errors bars in R p . The squares represent the SS99 data with error 
bars. 
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7. Tides 

The effect of tides was studied by several authors (Kantha et al., 
1995; Munk, 1966, 1997; Munk and Wunsch, 1998; St. Laurent 
et al., 2002; St. Laurent and Garrett, 2002; Garrett and Kunze, 
2007; Munk and Bills, 2007) and it requires the modeling of three 
distinct processes: (a) enhanced bottom diffusivity due to baroclinic 
tides, (b) tidally induced drag, and (c) unresolved bottom shear, 
which we now discuss in that order. 

7.1. Internal, baroclinic tides 

To generate bottom mixing, the key physical process has been 
identified to be the conversion of barotropic into baroclinic tides 
caused by the interaction of the former with rough bottom topog- 
raphy. The non-linear interactions among the baroclinic tides (and 
the shear they contain) allow part of their energy to be used to 
raise the center of gravity and thus produce mixing. 

The conversion of barotropic tides into baroclinic internal tides 
was studied by several authors, e.g., Kantha and Tierney (1997), 
Llewellyn Smith and Young (2002), St. Laurent and Garrett 
(2002) and Legg (2004) and was included in OGCMs by two groups 
(Simmons et al„ 2004, cited as S4; Saenko and Merryfieid, 2005, ci- 
ted as S5). Here, we employ the work of one of the present authors 
(S.R. Jayne). One begins by solving offline the 2D Laplace tidal 
equation with a resolution of 1/2° to obtain the barotropic tidal 
velocity u r , the 2D dynamical equations contain a drag which de- 
pends on the bottom topographic roughness denoted by h which 
is taken from the Smith and Sandwell (1997) data at 1/32° resolu- 
tion and then binned into the 1/2° resolution of the 2D code that 
provides u f . With the latter, one then constructs an expression 
for the internal tidal energy E(x,y) using the following parameter- 
ization by Jayne (2009): 

E{x,y) = \ pNKh 2 u 2 (W nr 2 ) (69) 

where (x',fi) are the wavenumber and amplitude. As discussed in 
Jayne (2009), the topographic roughness h 2 was derived from high 
resolution bathymetry [US Department of Commerce, 2006: 2-min 
Gridded Global Relief Data (ET0P02v2). National Oceanic and 
Atmospheric Administration, National Geophysical Data Center. 
Available online at http://www.ngdc.noaa.gov/mgg/fliers/06mgg01. 


html; Smith and Sandwell, 1997] as the root-mean-square of the 
topography over a 50-km smoothing radius, and k is a free param- 
eter set as k = 27t/125 km. It should be emphasized that Eq. (69) is a 
scale relation and not a precise specification of internal tide energy 
flux. In the barotropic tidal model, the value of k = 2rc/l 25 km was 
tuned to give the best fit to the observed tides. To construct the re- 
quired £ t ides. we employ the model suggested by St. Laurent et al. 
(2002) which has the following form: 

pfitides = qE{x,y)F(z) 

F(z)=Ar 1 exp[-(H + z)/a, 7T 1 = 1 - exp[-(H/Q] 

where the role of z~ ] is played by the scale function F(z) in which 
C = 500 m. The parameter q accounts for the fact that only a fraction 
q of the baroclinic energy goes into creating mixing; the remaining 
part 1 - q is radiated into the ocean interior where it may contrib- 
ute to the background diffusivity. The last step is the construction of 
the tidally-induced diffusivity using relation (lb): 

K a = r.(Ri,R p )^ (71) 

Since S4,5 also used (69)-(71 ), we need to point out the differences 
with their analysis. First, S4,5 used r = 0.2 for heat and salt while for 
momentum S4 took K rn = )0I< hs and S5 took K m = 10~ 4 m 2 s _1 , 
whereas in our model we have different heat, salt and momentum 
mixing efficiencies that depend on Ri and R p . This means that since 
these efficiencies are different, the mean T, S and velocity will be af- 
fected differently by tides. Second, we have updated the Jayne and 
St. Laurent’s (2001) original method. In particular, the model do- 
main was expanded to cover the global ocean (rather than ±72° as 
in the original work). Additionally, the gravitational self-attraction 
and loading in the tidal model was implemented in an iterative 
manner as in Arbic et al. (2004). Overall, these changes improve 
the fit of the simulated tides slightly: the diurnal tides improved 
significantly, likely due to including all of the Southern Ocean, 
where the diurnal tides are large, and the simulated semidiurnal 
tides did not improve. Though other parameterizations of the inter- 
nal wave conversion were suggested by Arbic et al. (2004) and Eg- 
bert et al. (2004), it was found that all of the schemes gave 
comparable accuracies in the simulated tidal elevations. The new 
expression for E(x,y ) is taken from Jayne (2009). Third, in the case 
of tides the Ri in (71) was taken to be the background value 0.5. 



Fig. 10a. Base 10 logarithm of the internal tidal energy flux E, Eq. (69), in units of W m 2 . 
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Fig. 10b. Drag power in W m 2 from the tidal velocities u, from Jayne and St. Laurent (2001). In Figs. 10a and 10b, the model is extended past ±72° that characterized the 
analogous figures in Jayne and St. Laurent (2001). 


7.2. Tidal drag, shallow seas 

Since the tidal energy of 1.51 Terawatts (Egbert and Ray, 2000, 
2003) dissipated as tidal drag is only 30% smaller than the one in 
internal tides, it is necessary to account for it. As is the case in 
Fig. 10a, the drag power in W m~ 2 shown in Fig. 10b corresponds 
to the updated model while the analogous figure in Jayne and St. 
Laurent (2001) corresponded to the old model ±72°. Contrary to 
internal tides, tidal drag cannot be represented by a diffusivity 
and its modeling is a non-trivial problem for several reasons. The 
bottom tidal velocity is generally larger than the mean velocity, 
for example, in shallow seas the tidal velocities are 0(5)cms~' 
which are much larger than 0(0.5) cm s _1 characterizing the mean 
velocities (Webb and Suginohara, 2001 a, b; Garrett and St. Laurent, 
2002 ). 

We begin by considering the component of the tidal field’s 
velocity that is along the direction of the mean field which can 
be modeled as an increased mean drag. That is done by extending 
the traditional quadratic bottom drag formula that depends only 
on the resolved mean flow u to include the tidal velocities u f so 
that the total velocity field is now u = u + u f . Beckmann (1998) 
and Haidvogel and Beckmann (1999, Eq. (5.19)) suggested the fol- 
lowing expression: 

t, = C d u|u|^C d u(u 2 + u2) 1/2 , C D = 0.003 (72) 

but since we did not find a derivation of it, we present one in 
Appendix C. The tidal velocities were taken from the same tidal 
model used to compute the function E(x,y ) in Eq. (69) and therefore 
the tidal contribution is location dependent. 

Two OGCMs have employed (72), the OCCAM Code 66 level 
model (SOC Inter. Report No. 99; http://www.noc.soton.ac.uk/jrd/ 
occam, 2005) and the GFDL Code (Griffies, private communication, 
2008). However, in both cases, the tidal velocity was taken to be 
constant while we employ the one derived from a tidal model 
and thus location and topography dependent. 

7.3. Unresolved bottom shear 

The component of the tidal field not aligned with the mean 
velocity cannot be modeled as a tidal drag. Since its mean shear 


is not zero and often large, it gives rise to a large unresolved shear 
Zunr with respect to the ocean’s bottom. This additional shear de- 
creases the local Ri possibly bringing it below Ri = 0(1), thus allow- 
ing shear instabilities to occur which ultimately enhance the 
diffusivities. 

We know of only one work (Lee et al., 2006) that includes T unr 
in an OGCM using a heuristic expression for Z^nr that depends on 
the M2 tidal velocity obtained from satellite data (Egbert et al., 
1994). Rather than using a heuristic expression for i unr , we 
adopted the viewpoint that since modeling an unresolved shear 
is a problem that has been widely studied in the context of the 
PBL (planetary boundary layer), there is a well assessed formalism 
we can adopt and which results in the following expression 
(Businger et al., 1971; Kaimal and Finnigan, 1994; Cheng et al., 
2002 ): 

^unr — ~~ (73) 

where u* is a frictional velocity and <P m (z/L) is a dimensionless 
structure function of the dimensionless ratio z/L where L is the 
Monin-Obukov length scale. Several field tested expressions for 
<T>{Ri) are available in the literature (e.g., Kaimal and Finnigan, 
1994). However, since Cheng et al. (2002) derived an expression 
for <P(Ri) from the RSM which is the formalism used in this work, 
for consistency reasons, we have adopted Cheng et al.’s expres- 
sion for <P(Ri) which was shown to reproduce previous empirical 
forms assessed against field experiments, the classical one being 
the Kansas experiment discussed in detail by Businger et al. 
(1971). As for u*. we model it in terms of the mean and tidal 
velocities. To do so, we consider the first relation (72) with 
u = u + u t : 

it,(total) = C D u(u 2 + 2u - u f + u 2 ) 1/2 

+ C D u t (u 2 + 2u ■ u t + u 2 ) 1/2 (74) 

In order to exhibit the contribution of the unresolved scales, we 
subtract from (74) its average thus yielding the unresolved part 
which then gives the desired u *: 

r f 1/2 

itfunr) = Tf,(total) - Tf,(total) u 2 t = Tjj(unr) (75) 
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In Appendix C we derive the following expression: 

[^(unr)] 1/2 = C d (u2) 1/2 (u 2 +u2)’ /2 


the type N 2 (z) = Nle (z ^ /h (Zang and Wunch, 2001), which gives, 
using relation (80): 

( 76) w* = h^K p Q_{Sv)=Aw* =Ah i K p (83) 


The u 2 from the tidal model averaged to the resolved scales charac- 
terizing the OGCM one employs, is used in (76) and the results are 
substituted into (75) and finally (73) to construct the Richardson 
number in which the shear is now given by: 


Ri = 


N 2 

I 2 ’ 


X 2 = Zt 


(77) 


where Z 2 es is the square of the shear of the resolved velocity field. 


8. Diapycnal velocity 


Once an OGCM is run using the mixing scheme just presented, 
the resulting large scale fields can be used to evaluate the diapyc- 
nal velocity w* which is an important part of the discussion on the 
origin of the MOC (meridional overturning circulation, Munk, 1966, 
1997; Munk and Wunsch, 1998, cited as MW; Doos and Webb, 
1994; Doos and Coward, 1997; Toggweiler and Samuels, 1998; 
Webb and Suginohara, 2001 a, b). 

Since our mixing scheme includes DD processes and since we 
were unable to find an expression of w* that includes different heat 
and salt diffusivities, we present such formula with some compar- 
ison with previous expression and some qualitative implications. 
Multiplying the mean temperature and salinity equations by a rs , 
respectively, and subtracting the two equations, one obtains the 
following expression for the diapycnal diffusivity w*: 


Arw* = g 


0 8T 
aT dz[ I<h dz 


8 8S 
~ as 8z\ s 8z 


4^N 2 )-N 2 (1-R P )-^ g-r 


(78) 


where a z = 6a/8z and where the spatial variation of the coefficients 
a TtS is due to the non-linearity of the seawater equation of state. In 
(78) we have not included cabbeling and thermobaricity which can 
be added, as shown in Klocker and McDougall (2010, in press). From 
the second form in (78), it is easy to check that if one takes: 

octs : z-independent (79) 

Eq. (78) reduces to the first term only which is the form of the 
advective-diffusive balance used by MW: 


where A is the surface of the ocean, z = 0 corresponds to the ocean’s 
bottom and z = H is the surface value and AT 2 3Af 2 /6z = li -1 . In the 
MW paper, the integral of N 2 (z) between 1 and 4 km was taken to 
be gAp/p=glCT 3 ; in our notation this corresponds to h = 1.3 km 
and thus we obtain: 


I< p = 1 cm 2 s- 1 : Q = 28 Sv, K p = 0.1 cm 2 s” 1 : Q. = 2.8 Sv (84) 


the first of which coincides with the case studied by MW. Next, we 
consider the contribution of tides. Using Eq. (70), the previous mod- 
el for N 2 , and the MW model, we obtain: 


w* = AT 2 


8I<pN 2 

8z 



(85) 


Depending on the relative sizes of the two scale heights, h, f, there 
may be an upwelling or a downwelling. For example, using 
£ = 500 m, as discussed previously, Eq. (85) implies that tides cause 
downwelling: 

Q(tides) < 0 (86) 


Next, we consider the effect of Double Diffusion. Using Eqs. (lb) and 
(27), the P = e relation given by Eq. (26) becomes: 

Ri-'r m -r p = l (87) 


Since a necessary condition for DD processes to exist is the absence 
of strong shear, we take Ri » 1 and thus r p = — 1 which means that 
the diffusivity becomes: 


K — r — — 
a p- j p n 2 ~ 



( 88 ) 


and thus: 


w* = AT 2 


8eF p 
8 Z 



(89) 


Depending on whether, in the DD dominated regime, e(z) increases 
or decreases with depth, we may have either upwelling or down- 
welling due to DD processes. The data in Figs. 13-16 of Kunze 
et al. (2006) do not allow us to draw a firm conclusion. 


9. Conclusions 


W = N- 2 ^(K p N 2 ) (80) 

Since MW further considered only positive K p > 0, it means that 
they did not include DD processes: thus, the MM model for w* does 
not include non-linearities in the seawater EOS nor does it include 
DD. On the other hand, if one consider the case: 

NoDD : K h =K s = I( p ^D 

Non-linearities in the seawater EOS : ct TS : z-dependent 
Eq. (78) reduces to 

NV.l(DNVD N ;(^-r->^) (82) 

which is Eq. (23) of Klocker and McDougall (2010, in press). 

Even without numerical computations, one can use the previ- 
ous relations to derive some interesting results concerning the ef- 
fects of DD and tides. Clearly, what follows has an illustrative value 
only. We first employ a constant diffusivity and a profile of N 2 of 


The complete mixing model which is composed of the RSM re- 
sults and different models for the dissipation e, is summarized in 
Appendix B. In the local limit, which is a justifiable approximation 
in a stably stratified regime, the RSM is fully algebraic and it only 
requires the solution of the cubic Eq. (41 ). The reason why it is pre- 
sented in a nested form is twofold: the physics of the various terms 
is easier to follow and from the numerical-computational view- 
point, nested relations are more advantageous. The physical aspect 
of the RSM is exhibited in relations (1) which show the key role 
played by the mixing efficiencies r a or by the structure functions 
S a which depend on Ri, R p and on the dynamical time scale which, 
in units of the mean shear, forms the dimensionless combination 
(rZ) 2 whose dependence on Ri, R p is obtained by solving the cubic 
equation (41 ) which yields (28). These functions S % and r a are dif- 
ferent for heat, salt and momentum, as shown in Figs. 1 and 2. We 
have assessed the validity of the RSM results based on produc- 
tion = dissipation by comparing predictions vs. measured data. 
The tests without an OGCM are as follows. 

Turbulent Prandtl number vs. Ri, Fig. 3c. There are abundant data 
that yield the ratio of the momentum to heat diffusivity vs. Ri 
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though only in the absence of DD processes. The data are repro- 
duced satisfactorily. 

Flux Richardson number, Rf vs. Ri, Fig. 3d. The RSM predictions 
reproduce the data without DD satisfactorily. 

DD processes. Past work by several authors showed how difficult 
it has been to construct a mixing model with DD + background tur- 
bulence. Our predictions provide a reasonable fit to the oceanic 
data, specifically: 

Mixing efficiency r h (Ri,R p ), Fig. 4. The data exhibit a clear depen- 
dence on Ri and R p and, as discussed in Section 5.8, previous mod- 
els based on linear analysis and laboratory data were not 
successful in constructing a DD model in a mildly turbulent back- 
ground. Laboratory data correspond to regimes with no shear that 
is, Ri -> oo, and their use in an ocean context is of doubtful validity. 

Bimodal e-distribution, Fig. 5a. The finding by SS99 of a bimodal 
kinetic energy dissipation e, namely that in the SF regime Ri > 1, e is 
an order of magnitude smaller than in the case of turbulence Ri < 1 , 
is reproduced rather closely. 

Heat to salt flux ratio r(Ri,R p ), Fig. 5b. In the regime of vanishing 
shear, the RSM predictions of the heat to salt flux ratio reproduce 
satisfactorily the data presented in Fig. 10a of SS99. 

NATRE mass diffusivity, Fig. 9 (St. Laurent and Schmitt, 1999). 
The model error range (obtained using data ranges forR p ), in most 
cases lies inside the data error range and in all cases intersects it. 

Tides. To describe the effect of tides, one must account for three 
distinct features: internal tides, tidal drag and the unresolved bot- 
tom shear. Internal tides were modeled in the same way as previ- 
ous authors via Eqs. (69)-(71 ) but the mixing efficiencies were not 
taken to be the same for heat, salt and momentum, rather, they 
were computed from within the model using relations (5) and 
the function E(x,y ) was taken from an updated model by one of 
the authors (Jayne, 2009). Tidal drag, which is most relevant in 
shallow seas, was only approximately accounted for in previous 
studies whereas we employ the results of the tidal model to com- 
pute the bottom drag rather than assuming a constant tidal veloc- 
ity, as done previously. As for the unresolved bottom shear, we 
have now included the tidal velocities not aligned with the mean 
velocities since they increase the shear, decrease Ri and enhance 
mixing. To model the unresolved shear, we adopted a procedure 
that has been successfully used in PBL studies. 

To assess the effect of the physical processes described above on 
the ocean’s global properties, one needs to employ an OGCM with a 
relatively high vertical resolution. For example, tidal drag which is 
expected to be the strongest in shallow seas, cannot be well repre- 
sented in OGCMs in which some shallow regions are converted to 
land or deepened due to the coarse horizontal gridding and 
requirements of numerical stability. Moreover, the OGCM treat- 
ment of deep regions using very thick layers near the bottom 
may not be able to resolve the bottom boundary layer so that the 
new effects, which are highly localized to the bottom, as opposed 
to the tidal energy which radiates upward with scale height ~1/ 
2 km, may not be allowed to act in full. Furthermore, the fact that 
often OGCM assign only one depth to each gridcell, whereas in re- 
gions of rough topography the true depth of the ocean bottom var- 
ies greatly over a gridcell’s area, degrades the performance of the 
tidal model. OGCMs with both high horizontal resolution and finer 
spacing in the vertical near the bottom and/or a parameterization 
which accounts for the actual distribution of bottom depths within 
each gridcell, are needed to fully assess the new mixing scheme. 

Future studies should also include into the mixing scheme the 
effect of mesoscales (30-100 km) and sub-mesoscales 0(1 km). It 
is well documented that they both re-stratify the mixed layer thus 
producing a reduction in the mixed layer depth (Oschlies, 2002; 
Mahadevan et al., 2010). No OGCM that we know of has included 
such mixed layer effects since there is still no satisfactory param- 
eterization of such processes. In addition, suggestions have been 


made that even in the deep ocean mesoscales may not move 
strictly along isopycnal surfaces, as assumed thus far: if they do 
not, there is a further contribution to the diapycnal diffusivity in 
addition to the small scale one discussed here. Recent studies (Tan- 
don and Garrett, 1996; Eden and Greatbatch, 2008a, b) have con- 
cluded that the effect may not be negligible especially at the 
ocean bottom. 


Note added in proofs 

It has been pointed out to us that in a recent paper, T. Decloedt 
and D.S. Luther (2010) have suggested a polynomial alternative to 
the exponential scale function F(z) in Eq. (70), which, as seen in 
their Fig. 4, provides a better fit to the data. Such a new formula, 
their Eq. (2), can easily be adopted in this model instead of the 
function F(z) in Eq. (70). We want to thank Dr. J. Richman for point- 
ing out this reference to us. 
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Appendix A. ID form of the Reynolds stress equations 


Vertical heat flux, / = wO: 


Df 

Dt 


= -w 2 T z +g(a. T 8 2 - ct s 6s^ - T 1 7t 4 1 / 
Vertical salt flux, J 5 = ws: 


= -w 2 S / + g(a. T 8s - a s s 2 ) - T 
Temperature variance, 0 2 : 


D(r 


= -2/"T z -2t- 1 ti 5 1 0 2 


Dt 

Salinity variance, s 2 : 

Ds 2 

W 


= -2/S, z - 2i 1 7T 3 1 s 2 


(A. 1 ) 


(A.2) 


(A.3) 


(A.4) 


Temperature-salinity correlation, Os: 

^ = -/S z -/r z -rV(A (A.5) 

Traceless Reynolds stress tensor by = Xy - 2<5y7<T/3 (i, j = 1,2,3): 

Dby 8K 1 7 1 5, 

Dt “ 15 Si;/ 2 Zij + 2 BiJ x bii (A6) 


with: 


Zy = b ik V jk + bj k V ik , By = g (yf + V? - \ SyAjJ^j (A.7) 

where Sy = 1 /2(U,j + Up) and Vy = 1 /2(U,j - Up) are the (mean) shear 
and vorticity tensors and Jf is defined as follows: 

if = <x T jt - Ojf, k = ~(gpr% (A.8) 
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Since Eqs. (A.6) involve also the horizontal heat and salinity fluxes 
via the buoyancy tensor By, one needs to account for their presence. 
The corresponding equations are: 

Horizontal heat flux,/? = ufi, i = 1,2: 

jjjr = -UiWT, z —J h 8 z Ui - 

Horizontal salt flux,/? = UjS, i = 1,2: 

^ = -ujws z -fd z Ui - x-' 7t, 7 ; 

Dissipation-relaxation time scales: 

For R p > 0 and Ri > 0, 

, 7I 4 = 71° ( 1 + 


q / 


Ri 


1 + qR, 


(A.9) 


(A.10) 


7I 2 = 7t 2 °(l+Ri)- , [l+2Rii? p (l+^)- 1 ], 715 = 71°, (A.ll) 

71? = 71° = (27/Co 3 /5)- 1/2 (1 + (rp 1 )- 1 , 

71? = 1 /3, 7I 3 = 71? = 71? = 0+ 

where a = 10, Ko = 1.66 and a t was taken from Eq. (10b). For R p ^ 0 
and Ri > 0, we further have the relations: 


7ll,4 = 7T° 4 (1 + Ri) \ 712,3,5 = 71? 35 

On the other hand, for Ri + 0. n k = 7t? for any k. 


(A.12) 


Appendix B. Complete mixing model 

Here, we summarize the complete form of the mixing model. 
Diffusivities (a = heat, salt, momentum, density) 


2 j( z p i 

General form: IQ = S x — = T z -= , r x = -Ri(xZ) 2 S x 

C ^ 

w 2 

Structure functions : S a = A y — 

Heat and salt: 

A h = 7i 4 [l + px + 7t 4 7r 2 x(l - r _1 )]~\ A s =A h (rR p )~ 1 

Heat-to-salt flux ratio: 

_ oc T w0 7t 4 1 1 + qx 
— 0C s WS 71, Rp 1 + px 

Momentum: 


5m — An, , A m — 


Ami 

A m 2 


where: 


Ami — - 


714-71, + ( 7t, - — )(1 -r ’) 


xA h 


A m2 = 10+ (7T 4 - 7t,i?p)x + — (tX) 


Ratio 'f\ 


I( 


1 +A X + 17 -^m{xlf 


, X={\-r~ 3 )xA h 


15 10 

Dimensionless variables x, p and q : 

X=R!(tX) 2 (1 -Rp)-\ p= 7I 4 7t 5 - 7I 4 7I 2 (1 + R p ), 
q = 71, 7I 2 (1 + R p ) - 71, 7l 3 R p 

Dynamical time scale C m = {xZ) 2 in the P= e, model: 
Cubic equation valid throughout the water column: 

C 3 G 3 + c 2 C 2 + c, G m + 1 =0 


(B.l) 

(B.2) 

(B.3) 

(B.4) 

(B.5) 

(B.6) 

(B-7) 

(B.8) 

(B.9) 

(B.10) 


with: 

c 3 = AjRi + A 2 Ri , c 2 = A 3 Ri +A 4 Ri, c, =A 5 Ri + Ag (B.l 1 ) 
where: 

150(1 - R P ) 3 A] = 7I,71 4 (7I 4 - 7i,R p )|7i 2 (157t 3 + 7) (r 2 + l) 

+ [14(7i 2 - 7t 3 ) - 157I?]R p } 

9000(1 - R P ) 2 A 2 = 7t,7i 4 {7t 2 (2107i, - 150 ti 3 + 1){r 2 + 1 ) 

+ [14(7i 2 - 7r 3 ) ( 1 + 1571, + 157t 4 ) + 15071 ?]R p 
+ 2107i 2 (7r 4 - 71 , )| 

150(1 - R p ) 2 A 3 = 7i, [57i 2 7t 4 (307t 3 + 17) + 7r,(157t 3 + 1)}{r 2 + l) 

- (157T 3 + 7)(7I? - 7 I 4 ) - [1 07T, 7t 3 7T 4 (1 57T 3 + 17) 
+ 157t 2 (7t? + 7l|) + 147T, 7T 4 (1 - 107I 2 )]R p 

9000(1 - R P )A 4 = [150(7t, 7t 3 + 7i 2 7i 4 ) - 771,(1 + 307r,)]R p 

- 1 50(71, 7I 2 + 7t 3 7I 4 ) + 7 tt 4 (1 + 307I 4 ) 

30(1 - R p )A 5 = [-30(71, 7t 3 + 7t 2 7i 4 ) - 1771 , }Rp 
+ 30(71, 7I 2 + 7t 3 7t 4 ) + 17714, 

Ag = -l/60 (B.12) 

The n’s are given by Eq. (A.ll). 

Dissipation 
Mixed layer: 

e = B 2 (xir 3 £ 2 Z 3 , t = -jji-J , 

4 = Kze 0 (£ 0 + kz)~' (B.13— B.15) 

where ( 0 = 0.1 7H, k = 0.4 is the von Karman constant and H is the 
ML depth computed as the point where the potential density and 
the surface value differ by jcr(H) - er(0)| > 3 x 10~ 5 g citT 3 ; R f is de- 
fined in Eq. (42) and B, = 21.6. 

Thermocline: 

£ = £ igw L(0,JV) (B.l 6) 

the dimensionless function L(0,N) is given by Eq. (65b), the expres- 
sion for gjgw N~ 2 is taken from the Gregg-Henyey-Polzin model, Eq. 
(66). 

Tides 

pfitides = QE(x,y)F(z), E(x,y) = ^ pNKh 2 u 2 (B.17) 

where the wavenumber k = 2tc/ 1 0 km and h is the roughness scale 
obtained from Smith and Sandwell (1997). The scale function F(z ) 
has an exponential shape with a spatial decay scale of £ = 500 m: 

F(z)=AQ 1 exp[-(H + z)/{], A-’ = 1 - exp [-H/Q (B.18) 

Bottom drag 

x„ = C D u(u 2 + u 2 V /2 , C D = 0.003 (B.19) 

Unresolved bottom shear 

In the BBL, the Ri that appears in Eq. (B.l ) must be taken to be: 

= X2 = 4s + X u„r> X unr=^<Z>m, 

/ \ 1/2 / „ \ 1/2 

uJ = C D (u 2 J (u 2 +u 2 J (B.20) 

where the function <P m is given by Eq. (36) of Cheng et al. (2002). 
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Fig. 11a. Area weighted histogram of the ratio of the rhs to the ihs of the relation for the mean magnitude of the unresolved stress, Eq. (C.8). For details, see text. 



Fig. lib. Same as in Fig. 11a but the histogram is weighted by the product of the area and the lhs of Eq. (C.8). 


Appendix C. Relations (72) and (76) 

The total velocity field is contributed by mean and tidal 
velocities: 

i),(total) = C D u|u|, u = u + u t (C.l) 

We use an overbar to denote the averages used in OGCMs, 
u f = 0 and u t |u + u t | = 0 since by symmetry the latter vector can 
only point along the direction of u and, to the extent that the mean 
and tidal fields are uncorrelated (because they represent low and 
high frequency fields), we expect such a term to vanish. We then 
obtain: 

T ft (total) = C D u(u 2 + 2u ■ u t + u 2 ) 1/2 (C.2) 


Next, we neglect the second term in the parenthesis since it is the 
product of high and low frequency variables with little overlap thus 
giving a zero mean. Eq. (C.2) then becomes: 

T^total) = C D u(u 2 + u 2 ) 1/2 (C.3) 

If one exchanges the square root with the averaging process, one 
obtains Eq. (72). Numerical experiments by Saunders (1977) sug- 
gest that even in the worst case the error in making this approxima- 
tion is no more than 50%. By contrast ignoring the tidal contribution 
to the drag altogether, as done in most OGCMs to date, can lead to 
an error of an order of magnitude. As for Eq. (76), we begin with Eq. 
(C.l ). Writing u 2 for ii u, we have: 

t b (total) = C D (u + u,)(u 2 + 2u ■ u t + u 2 ) 1/2 


(C.4) 
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In order to exhibit the contribution of the unresolved scales, we 
subtract from (C.4) its average yielding the unresolved part: 

Tf,(unr) = T;,(total) - T b (total) (C.5) 

from which we derive u* as follows: 


uj = [T2(unr)] 1/2 (C.6) 

Thus, the procedure consists of taking the modulus of (C.5) and 
averaging it over the OGCM scale. Even adopting the approxima- 
tions used in (C.3), the expression for T^(unr) turns out to be still 
rather complex: 

C D 2 T^(unr) = (u 2 + u 2 ) 2 + u 2 d 2 - 2u 2 (u 2 + u 2 ) 1/2 d, 

A = (u 2 + u 2 ) 1/2 (C.7) 

A similar approximation to the one that led from (C.2) to (70), 
yields: 

u 2 t = jr 2 (unr)] ’ /2 = C D (uf + u 2 u?) 1/2 

=>C D u? 1/2 (u 2 +u 2 ) 7 (C.8) 

where the last approximation was necessary because we lack data 
on nf. The last relation is (76). 

Lacking access to fine time and space scales ocean velocities 
field data including mean and tidal components, we tested the 
above approximation using simulated data. The latter were created 
by interpolating the velocities from our 3x3° NCAR OGCM similar 
to that used in Canuto et al. (2004b) to the 1/2 x 1/2° grid of the 
tidal model and then adding at each tidal gridbox a linearly polar- 
ized sinusoidal in time velocity field with rms magnitude equal to 
that of the time-averaged tidal velocity square of the tidal model’s 
output. Four polarizations, east, northeast, north and northwest 
were used and 12 time steps were taken. We then computed the 
left- and right-hand sides of (C.8) from the simulated data for each 
of the polarizations, where the overbar was taken to be an average 
over time and the 3x3° gridcell. The rhs of (C.4) was substituted 
into (C.5) to compute T 6 (unr) which was then substituted into 
the lhs of (C.7). Finally, the ratio of the rhs to the lhs of (C.8) was 
computed for each polarization for each gridbox. The average 
weighted either by gridbox area or gridbox area times the lhs of 
(C.8) was 0.99. Histograms of this ratio with the former and latter 
weighting are presented in Fig. 11, respectively. We consider the 
results adequate empirical evidence of the validity of (76) in the 
context in which we are applying it. 
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